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ABSTRACT 

The process of diffusive shock acceleration rches on the efficacy with which hydromag- 
netic waves can scatter charged particles in the precursor of a shock. The growth of 
self-generated waves is driven by both resonant and non-resonant processes. We per- 
form high-resolution magnetohydrodynamie simulations of the non-resonant cosmic- 
ray driven instability, in which the unstable waves are excited beyond the linear regime. 
In a snapshot of the resultant field, particle transport simulations are carried out. The 
use of a static snapshot of the field is reasonable given that the Larmor period for 
particles is typically very short relative to the instability growth time. The diffusion 
rate is found to be close to, or below, the Bohm limit for a range of energies. This pro- 
vides the first explicit demonstration that self-excited turbulence reduces the diffusion 
coefficient and has important implications for cosmic ray transport and acceleration 
in supernova remnants. 
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1 INTRODUCTION 

Supernova remnants (SNR) have long been identified as pos- 
sible locations for the production and acceleration of galac- 
tic cosmic rays. The diffusive shock acceleration mechanism 
provides a natural explanation for the observed power-law 
spectrum for these cosmic rays. Acceleration to high en- 
ergies relies on confinement of relativistic particles to the 
accelerating region close to the shock. Pitch-angle scatter- 
ing from self-produced hydromagnetic waves can provide a 
suitable mechanism, but a detailed analysis demonstrates 
that the sp ectrum still falls short of the cosmic ray knee at 
lO^'^ eV ijLagage fc Cesarskvl[l983l '). However, theory dic- 
tates that the maximum energy of the cosmic rays is deter- 
mined by the diffusion coefficient, which is assumed to scale 
with the strength of the ambient magnetic field. Therefore, 
in the quasilinear framework, amplification of the apparent 
ambient magnetic field in the vicinity of the shock facilitates 
acceleration to higher energies. 

To within an order of magnitude the timescale for the 
diffusive shock acceleration of a particle is tacc ~ f^/Vsh 
where k = Iimt_oo{A^r-)/2t is the asymptotic spatial dif- 
fusion coefficient where Ar is the spatial displacement over 
time interval t. The conventional approximation for trans- 
port of particles of speed v via scattering is Bohm diffusion 
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which assumes a single scattering event occurs during each 
Larmor time. 

In the region upstream of a shock, the scattering of 
particles is due to resonant collisions with slowly mov- 
ing magnetic irregularities which, in turn, are amplified 
by such interactions themselves. These turbulent irregulari- 
ties, until recently, have been assum ed to saturate at levels 
SB < Bo (jMcKenzie fc VolklflQsi ). With such turbulence, 
Bohm diffusion can be taken as a lower limit for the diffu- 
sion coefficient, K>KBohm. In highly turbulent fields, how- 
ever, numerical investigations suggest that particle trans- 
port pr operties may differ significantly f rom this approxi- 
mation ijCasse. Lemoine fc Pelletierll20oj ). 

iBell fc Lucekl l)200l') proposed a process of resonant am- 
plification of the magnetic field driven by the pressure gra- 
dient of cosmic rays, allowing rapid transfer of energy into 
the Alfven waves. Amplification factors of SB /Bo ~ 1000 
were estimated from the linear theory, although this value is 
possibly excessive as the analysis most likely does not hold 
beyond SB /Bo ~ 1. Nevertheless, a more detailed non-linear 
analysis does suggest acceleration up to i?max ^ 10^'' eV. 
In particular, under expansion into a stellar wind, it may 
even be possible to achieve proton energies of the order 

-Em 



10 eV via this process. 



Moving beyond this result. [Belli ()2004l ) identified a non- 
resonant instability driven by streaming, positively charged 
cosmic rays that operates on length scales much shorter than 
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the Larmor radius rg = p/eBo of the driving particles (where 
p is the particle momentum, e is the electronic charge, and 
Bo is the mean magnetic field magnitude). Crucially, under 
typical conditions in the precursor of a high Alfven-Mach- 
number supernova remnant shock, the growth rate for this 
non-resonant instability can be several orders of magnitude 
greater than that of its resonant counterpart. Indeed, numer- 
ical magnetohydrodynamic simulations carried out in the 
same work confirmed that field amplification by a factor of 
an order of magnitude is possible. 

Only recently have analytical approaches to particle 
transport achieved results consistent with numerical ap- 
proaches. In the case of interplanetary shocks it has been 
shown that using detailed models of particle transport in 
tur bulent environment s sub-Bohm diffusion can be achieved 
fsee lZank et al.l l|2006l ) and references therein). 

In this paper we investigate, using high resolution MHD 
simulations, the non-resonant Bell-type instability and its 
non-linear evolution. We also determine the effective trans- 
port properties in the amplified field and compare them with 
the standard Bohm diffusion description. 

The structure of the paper is as follows. In Section[2] we 
recall the linear MHD analysis of the dispersion relation rele- 
vant to the numerical simulations performed. In Section[3]we 
describe the numerical method used to investigate the non- 
resonant instability and discuss some properties of the field 
produced. In Section 3] we describe the techniques used to 
integrate equations of motion in the snapshot field obtained 
from the simulations. We discuss the transport properties of 
both field lines and test particles. Finally, in Section O we 
present our conclusions. 



where s is an index used to sum over the charged species and 
/^o is the plasma permeability. The zeroth order cosmic-ray 
current density is parallel to the mean magnetic field, which 
is chosen to lie along the z-axis. All first order perturbations 
are taken to lie perpendicular to the zeroth order compo- 
nents such that the magnetic field and cosmic-ray current 
density can be expressed as B = By -|- Bx, jcr = jy + j± 
respectively. Making use of the induction equation 



f = Vx(uxB) 



(3) 



we look for plane-wave solutions to the linearised system 
with k also parallel to the mean field. All first order quan- 
tities Bx, jx, ux and Ex = — ux x By, take the circularly 
polarised form ± iAy = yloexp[j(fcz — (jJt)]- For the pur- 
poses of analysis, the plasma is taken to be initially at rest. 

Using a kinetic description, it was shown by 
iReville et all \2m% that for modes with wavenumbers 
kr^i 3> 1, where rgi = p\/e\B\\ \ is the Larmor radius of the 
cosmic rays of minimum momentum pi , kinetic effects can be 
neglected i.e. both thermal effects and perturbations to the 
mean cosmic-ray fiux are negligible. Splitting the perpen- 
dicular magnetic field perturbation into its cartesian com- 
ponents in the x-y plane the relation 



-iB\\3\\k/p 



up- - vlk^ 



By 



(4) 



is obtained where — B^^/y/jiop is the Alfven speed in the 
unperturbed field. This system permits several waves, whose 
dispersion relation takes the form 



2 THE NON-RESONANT INSTABILITY 



Following iBelll l|2004 ). we determine the dispersion relation 
for waves in the precursor of a quasi-parallel non-relativistic 
shock, where the energy density of the cosmic rays is compa- 
rable to the ram pressure of the shock (~ 10%). The usual 
assumptions made in the linear analysis of MHD waves prop- 
agating in a plasma are that the plasma is quasi-neutral, i.e. 
is charge neutral on aggregate. Therefore, assuming a plasma 
consisting entirely of protons and electrons in the presence 
of a pure proton cosmic ray component, the thermal plasma 
has a charge excess of electrons. In the reference frame in 
which the upstream protons are at rest, the cosmic rays are 
seen to stream with a group velocity approximately equal 
to the shock velocity, (v) ~ Vsh- A charge flux is induced in 
thermal plasma to neutralise this current. In the case of a 
mean field By parallel to the shock normal, this return cur- 
rent has a mean (unperturbed) value given by jy = ricrevsh 
where the number density of cosmic rays and e is the 

electronic charge. The momentum equation is given in the 
MHD approximation by 



du „ „ . „ 

p-^ = -VP + jthXB, 



(1) 



where jth, the current density carried by the thermal plasma, 
is determined by Ampere's law. 



V X B = fj.o'^nsqs-Vs = Mo(jth + jcr). 



(2) 



= k v^± (Vsh 



(5) 



similarly to equation (15) from lBell l|2004 ) where the dimen- 
sionless driving parameter ( is given by 



C = MO 



Bh 



-rgi. 



(6) 



It follows that, for strongly driven modes such that 
Cvsh/i^a 3> 1, there exist aperiodic waves (Re co — 0), over 
the range 1 < kr^i < C'^'sh/^a- From equation ([5]) it is 
straight forward to show that the fastest growing mode oc- 
curs at 



I, _ MO Jll 
with growth rate 



(7) 



(8) 



The dispersion relation for arb itrary orien tations of k, 
jcr, and Bo has been determined byHel l|2005h demonstrat- 
ing that unstable modes exist for all orientations provided 
k ■ Bo 7^ 0. Melrose (2005) has compared the resonant and 
non-resonant instabilities for arbitrary angle of propagation 
and shown that the latter occur with elliptically polarised 
modes. This emphasises the need for high resolution simu- 
lations in three dimensions. 
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Figure 1. The root-mean-square value of the perpendicular 
(solid) and parallel (dashed) components of the magnetic field. 
The field is initially damped as the initial energy input goes into 
kinetic energy, setting the inertial plasma into motion. Once the 
fluid is in motion the perturbations grow according to the linear 
analysis as exp(7inaxt). 



3 INSTABILITY DRIVEN TURBULENT 
FIELD GENERATION 

Our simulation box represents a region in the restframe of 
the upstream magnetised fluid (shock precursor). We assume 
a thin box such that both the pressure and density can be 
taken to be initially constant, and the cosmic rays can be 
considered unmagnetised (i.e. they travel through the box 
with straight trajectories ). 

Following iBeil (|2004l ) we take typical interstellar values 
for the upstream quantities By = 30nT, nucleon density 
no = 1 X 10~^m~'^, corresponding to an Alfven speed Wa = 
6.6x lO^ms"^. Furthermore assuming 10% energy transfer to 
cosmic rays and shock speed Vs/c — 1/30, with ln(p2/m-c) = 
14, where p2 is the upper cut-off of the cosmic-ray energy 
spectrum, we obtain ^ — 4.8 x 10"** and (^v^ — llOOua. 

In the quasilinear theory of diffusive shock acceleration, 
the minimum momentum of the cosmic-ray distribution is 
expected to increase with distance upstream, since lower en- 
ergy cosmic rays are confined closer to the shock (Eichler 
1 19791 : lBlasill2002l ). Assuming Bohm diffusion the minimum 
cosmic ray momentum at a given distance z upstream of the 
shock is pi ~ SVsheBz/c. In order for the non-resonant mode 
to leave the regime of linear growth before being overtaken 
by the shock front, i.e., before being advected over a distance 
of roughly crgi/vgh, we take pic = 1 PeV (although this is 
not a lower limit). This gives corresponding physical scales 



k^L = 2 X 10^=* m. 



Tmax 



= 98 yr, and r^i = 1.1 x 10' 



3.1 Turbulent seed field 

The non-resonant instability develops in a field containing 
a level of seed turbulence. In the simulations that follow 
the magnetic field is initialised according to the following 
prescription 



B(a;,2/,2) = Bo +SB{x,y,z), 



(9) 



where 5B is a small isotropic turbulent field component 
{SB/Bo <C 1). The turbulent c ompo nent is constructed sim- 
ilarly to iGiacalone fc Jokipiil l| 19941 ) , by choosing A'^ plane- 



wave modes with random phases, orientations, and polar- 
ization states given by 



5B{r) = 2^A„ 



(10) 



where kn = fcncj,. For each mode n, A„, ipn, k„, and ej, 
are the amplitude, phase, wave-number, and propagation 
direction respectively. Additionally, the polarization vector 
is given by 



^„ = cos(a„)e„ + i sin(Q:„ 



(11) 



where Qf„ is the polarization angle. The three vectors (e^, 
e^j, efj) form an orthonormal basis for the plane waves such 
that, for a continuous representation, the field is guaranteed 
to be solenoidal by the condition kn -^n ~ 0. However, when 
the field is projected onto a mesh of uniform spacing h the 
solenoidal condition on a cen tral-differenced grid is given by 
lO'Sullivan fc Downed (|2007l ). 



5n • sin(k„/i) = 0. 



(12) 



The amplitude of each mode for a given variance is 



Ai^ = cr^Gn 



Eg. 



with 



Gn — 



AV„ 



l + (fcnic)r 



(13) 



(14) 



Here F is the spectral index of the turbulence, and the nor- 
malisation factor for three-dimensional turbulence is AV„ — 
4-Kk'^Akn- Lc is the correlation length of the magnetic field, 
which in this case we take to be the largest wavelength in 
the system. 

In order to generate the turbulent field component 4 
random numbers are required for each mode: two to specify 
the orientation of the wave- vector e^, and one each for the 
phase and polarization angles tpn and Qn respectively. 

Since the turbulent field will be represented on a peri- 
odic mesh, each mode must have an integral number of wave- 
lengths in each coordinate direction. Explicitly, assuming a 
cubic grid of dimension L with A'' mesh points in each coor- 
dinate direction, k = 2TT/L{nx,ny,nz) where rix = 0..N/2 
etc. (The factor 2 is necessitated by the hermiticity of (5B 
requiring both negative and positive k modes be allowed.) 
This geometrical constraint also means that at small k, the 
density of unique modes is low and available k vectors will 
have a strong preferential bias along the grid axes. Since 
isotropic and homoge neous turbulence is only achieved in 
the limit of A'^ — » oo (|Batcheloij[l98l ). we do not include 
modes from the low-fc range where this bias is greatest. Ad- 
ditionally, we attenuate the dynamic range at high-fc by re- 
quiring the shortest scale fluctuations are well resolved over 
10 cells (rather than 2 as minimally required above). 



3.2 Numerical Scheme 

We describe the numerical scheme used to perform the simu- 
lations of the non-resonant Bell-type instability. The simula- 
tions are performed with consid erably high er resolution than 
those performed previously bv iBelll l|2004h . The cosmic-ray 
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current is taken to be uniform along the 2-axis and con- 
stant in time. Kinetic effects are again neglected such that 
jx = 0. For clarity from this point onwards, terms involv- 
ing magnetic field B are assumed to have absorbed a factor 
l/y/7o, and terms involving charge current density j are as- 
sumed to have absorbed a factor ^/lo. The MHD equations 
including cosmic ray current may be written in the form 



9U 



-HV- V 



dt 
where 

U = (p,pu,B,£)^ 
and 



(15) 



(16) 



V = 



(17) 



( 

puu + Ip* - BE 
uB - Bu 
V (£;+p*)u- (u-B)B / 

with p* = p + the total pressure. The source terms 

are given by 



S = 



^ ° ^ 

-j||Z X B 



V -u- (ji|Z X B) / 



(18) 



We have constructed a finite volume Godunov code in 
three dimen sions, MENDOZA , which follows the method 
described in iFaUe et all (Il998l ). The equations are split into 
three one-dimensional problems to find the fiuxes on the 
surface of each cell. The fiuxes are determined at each cell 
interface from the solution to the linearised Riemann prob- 
lem, and the field is updated with second-order accuracy in 
both time and space. The source terms are then added as 
explicit volume- averaged quantities. We have investigated a 
variety of divergence cleaning methods, but found t he gener- 
alised Lagrangian multiplier (GLM) method of D cdner et al.l 
l|2002l ) to be most effective. Divergence errors are naturally 
produced in multidimensional simulations, due to the fact 
that the one-dimensional problem is solved for each direc- 
tion and the fluxes simply added together. While the GLM 
method does not identically conserve V • B = 0, it has the 
unique property of rapidly damping any divergence intro- 
duced, while advecting them away from the problem areas 
at the fastest admissible speed. For the study of magnetic 
instabilities, particularly in the presence of strong pressure 
gradients, we find this method to have many advantages 
over correction schemes. For a comparison of many differen t 
divergence cleaning methods see, for example, iTothI l|2000l ) . 

Using simulation units, the model is then initialised 

with 

P = p=l, j|| = 47r/5 , B|| = 1 , u = . 
Equations ((SI to ((8| then yield 

t'a = 1 , «sh = 1520 , fcmax = 2-7r/5 , 

7max = 27r/5 , rgi = 1375/7r . 

As is common practice in the literature, the turbulent 
field spectrum is represented by a superposition of plane 
wave modes, one selected in a random direction from each 
of a number of uniform intervals of log k over a finite range. 
As commented in S ect ion [3 . 1 1 however . since a discrete mesh 
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Figure 2. Slice through the x-y plane at the centre of the simula- 
tion box at i = 6 showing the isoUnes of the magnetic field mag- 
nitude. Contours have uniform spacing of 0.5 with Smin = 0.5 
and Bmax = 13. The regions of strongest field are found in the 
regions of higher density. 



is being used to represent a periodic field, viable fc-vectors 
are finite in number and may be clustered or multiply degen- 
erate in k. Therefore, for a given grid resolution, requiring 
that each interval contains at least one viable fe-vector de- 
termines the total number of modes used to construct the 
turbulent field representation and the cutoffs in the power 
spectrum. In this work, a 512'^ grid is used with modes in 
the range 40 ^ kL ^ 550. With these parameters the num- 
ber of modes we obtain \s N = 151. While this dynamic 
range (~ 14) is relatively small, it does represent isotropic, 
homogeneous turbulence well over the chosen narrow band 
since at least one mode exists within each interval (by con- 
struction). As long as fcmax is central to this range, it should 
make little difference to the experiment that the spectrum is 
narrow since for k > 2fcniax the linear growth rate is zer o and 
for k < 2A;max the growth rate goes as fc^/^ (|Bellll2004l '). We 
choose L — 64 such that feniaxi ~ 80 and the corresponding 
fastest-growing fiuctuation scale is well resolved over ~ 40 
cells. 

In simulation units, the initial seed field is set with — 
O.OlBji. Finally, we restrict our studies to the spectral index 
corresponding to three-dimensional Kolmogorov turbulence, 
r = 11/3, (e.g. ,Giacalone fc JokiDii.1999, '). 

3.3 Results and Discussion I 

Initialising the fields using the parameters described in sec- 
tions [3]T] and |3]2l the development of the root mean squared 
parallel and perpendicular components of the magnetic field 
is plotted in Figure [T] Since the plasma starts from rest, the 
energy input at the beginning primarily goes into kinetic 
energy of the background fluid and no field growth is ob- 
served. After this brief initial stage, at approximately t — 1, 
the magnetic field starts to grow at close to the linear growth 
rate 7max. The energy in the perpendicular magnetic field 
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Figure 3. Logarithmic plot of tlie power spectrum in the am- 
plified field at times t = 5 and t = 6. The spectrum at large k is 
a power law with P{k) = fc"", with s x 7.5 at t = 5 and s ~ 5 
at t = 6. The turnover in the spectrum moves to smaller k as the 
parallel component of the field increases. 



eventually overtakes that of the parallel field but continues 
to grow at the same rate. As the amplitude of the unsta- 
ble modes approaches that of the mean field, cavity s truc- 
tures similar to those observed in the simulations of iBelll 
l|2004 ) start to appear. The development of these cavities is 
driven by the Lorentz force exerted on the circular/elliptical 
modes by the return current. The expansion of each cavity is 
hindered by the growth of neighbouring cavities and dense 
regions are formed between them. The frozen-in magnetic 
field becomes very large inside these cavity walls, as shown 
in Figure [2] 

It is at this point that neglecting the non-linear feed- 
back of the field o n the cosmic rays becomes apparent. As 
discussed in [Bdil (|2005l ). information is only passed along 
the direction of the cosmic ray current via tension in the 
magnetic field lines. However, as the field increases, the cos- 
mic rays themselves would most likely be beamed into the 
resulting cavities. From the simulations, it is observed that 
there are no clear coherent structures forming along the di- 
rection of the cosmic ray current. 

A second difficulty is that vacuum states begin to form 
within the cavities, and the MHD equations lose their valid- 
ity. A possible solution to this is to use a highly dissipative 
scheme, but this reduces the accuracy of the solution. Using 
such a method the results obtained for the l o ng ter m evolu- 
tion of the field are similar to those in iBelll l|2004 ). and we 
do not comment on them. 

Since we are primarily concerned with determining dif- 
fusion statistics in an amplified field we choose to take a 
snapshot of the field in the early stages of the non-linear 
development, at t = 6. The power spectrum of the field at 
t = 5 and t = 6 is plotted in Figure |3] 

In the shock frame, the amplified field is advected to- 
wards the shock. It is expected that lower energy particles 
are confined closer to the shock by the turbulence produced 
further upstream by the cosmic rays which can diffuse fur- 
ther from the shock. We determine the transport properties 
for relativistic protons in a static amplified field taken at 
t = 6. 




Figure 4. A comparison of the fields after being passed through 
the low power filter. The same large scale structure maintained 
quite well between the filtered field (left) and MENDOZA's field 
(right). 



4 COSMIC RAY AND FIELD LINE 

TRANSPORT IN TURBULENT FIELD 

We now discuss the first results of cosmic-ray transport in 
a turb ulent field resulting from the non-resonant instabil- 
ity of Hel l|2004l ) as described in Section [2] In particu- 
lar, we investigate how the transport properties differ from 
the Bohm diffusion limit. Much work has been done on 
this topic numerically by llGiacalone fc Jokipiilll99"3 , 1 19991 : 



ICasse. Lemoine fc Pelletierl 20021 . for example) but in all 
cases an a priori field was assumed. However, the field we 
use has been determined self-consistently. 

To determine the statistical transport properties 
of the amplified field we h ave developed the code 
"LYRA" (|0'Sullivan et al.ll2008h . LYRA assumes a Fourier 
mode description of a static magnetic field as described in 
Section [3.11 However, MENDOZA represents the magnetic 
field at a finite number of mesh points (in our case, at 512"^ 
uniformly spaced points). We convert the field into a contin- 
uous Fourier mode representation by taking a Fourier trans- 
form and importing the 5000 highest power real modes into 
LYRA. In the particular case considered here, passing the 
field through what amounts to a low-power filter results in 
a loss of 16% of the power in the field. Any lost power is 
returned to the field by rescaling the included modes ac- 
cordingly. Figure|4]provides a comparative illustration of the 
pre-processed field generated by MENDOZA and the filtered 
field read in by LYRA for the particle transport simulations. 
As can be seen, the overall large scale structure of the field 
is maintained quite well, but there is a certain amount of 
smoothing of the sharper features. However, the continuous 
description of the field provided by the Fourier transform 
is more practical than a piecewise constant mesh for the 
purposes of integrating particle trajectories. The resulting 
field is also divergence free. We can not be certain that the 
spreading of the field inside the cavity walls does not affect 
the final results, but in the current work the Fourier descrip- 
tion makes the computations much simpler. We will address 
this issue in future work. 

Another important limitation introduced by importing 
the field from a periodic grid representation is that while, 
in theory, we have access to an arbitrarily large dynamic 
range, the field will remain periodic on the relatively small 
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Figure 5. The Perpendicular magnetic diffusivity plotted 
against distance along a field line. Both in in units of box lengths. 



scale of the parent field. It is important to consider this point 
when interpreting results displaying any similarity on scales 
corresponding to the periodicity of the field. 

With the field described in this manner, the equa- 
tions of motion are integ rated using the Burlisch-Stoer 
method (jPress et al.lll993 ). In the case of field lines, the 
equation to be integrated is dr/ds = B/B. For protons, the 
governing equations are 



ev X B, 



(19) 



(20) 



d7v 
"dT 
and 

dr _ 

dt ^ 

When integrating the equations of motion for protons, the 
tolerances in the integrator are set to conserve energy to 
< 1% over 1000/27rtJ~^ where = 27rrg. We determine 
the statistical transport properties for 9 different energy val- 
ues ranging from 7 = 10"^ to 7 = 4 x 10^ . For each value, an 
ensemble of 2000 particles is initialised with random starting 
positions and pitch angles. Using the t = & field snapshot 
described in the previous section, with magnetic field simu- 
lation units expressed in terms of a 3/iG ISM field value, we 
carry out particle transport simulations using LYRA. 



4.1 Results and Discussion II 

To investigate the properties of the field lines themselves 
2000 field lines are integrated, and the mean square perpen- 
dicular displacement along each field line is calculated. The 
magnetic diffusivity D = {Ax'^)/2s, where s is the distance 
along a field line, is plotted in figure [5] Since the fastest 
growing modes are circularly polarised, the field lines are 
helical, although the guiding centre wanders about the x-y 
plane as it moves along the z-axis. The parallel diffusivity is 
essentially ballistic (Dy oc s). 

We numerically determine the statistical average diam- 
eter of the helices by iterating over the same 2000 field lines, 
and approximating the diameter of each revolution. In what 
follows, all length units are normalised to simulation box 
lengths. The root mean square diameter of the helices is 
d = 0.148 which is an important length scale in the study 
of the fields properties. The peak in the perpendicular mag- 
netic diffusivity is located at s = 0.273. Using this value 
and setting the perpendicular displacement equal to the root 




Figure 6. 



at t=6 showing times for super-difl'usive regime 



in units of box-lengths squared per Larmor time. Lines correspond 
to gamma factors (from bottom to top at At = 100) 7 = 4 X 
10^, 10"*, 2 X 10^, 2 X 10"', 10^. Time units are Larmor times 
for a 7 = 10* particle. Note these lines continue to asymptotic 
values and we show only the early evolution to demonstrate the 
super-difl^usive behaviour. 



mean square diameter gives D± =0.08 which is in reason- 
able agreement with the experimentally determined peak at 
approximately D± = 0.057, as shown in figure [S] 

The mode with the fastest growth rate has a wavelength 
similar to the gyroradius of protons with a particle gamma 
factor of 7 « 2 X 10'*. We focus on the transport properties 
close to this value. Particles with 7 much less or greater than 
this have mean free paths close to a cell width or box length 
respectively. 

We interpret these results in terms of the diffusion coef- 
ficients for the actual particle trajectories determined from 
LYRA. The parallel and perpendicular diffusion coefficients 
are determined by averaging over a large ensemble of parti- 
cles 



Kj_ 



lim 



lim 



2At 

(Axj) 
2At 



(21) 



(22) 



Note that the and k± refer to parallel and perpendicular 
directions with respect to the initial field, whereas in the 
t = 6 snapshot of the field the perpendicular component 
already dominates, as is evident from Figure [1] 

The initial evolution of the particle trajectories are gov- 
erned by the field line statistics. Figure [6] and Figure [7] 
plot the early parallel and perpendicular diffusion respec- 
tively, showing an initially super diffusive regime, with a 
peak in the perpendicular value for the lower energy parti- 
cles occurring at approximately 5 itJg^ ■ Inserting this values 
into Eq. (|22p and equating to the experimental peak value 
~ 0.0015, gives an estimate for the root mean square per- 
pendicular displacement at the peak, ^ {Ax'^) = 0.08. This 
value is also in agreement with the statistical mean radius of 
the helices. For higher energy particles the Larmor radius of 
a particle becomes comparable to or larger than the average 
diameter, and the peak time differs from this value. These 
particles are not trapped to individual field lines and must 
travel a greater distance before becoming diffusive. 
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Figure 7. k± at t = 6 showing times for supcr-difTusive regime in 
units of box-lengths squared per Larmor time. Lines correspond 
to gamma factors (from bottom to top) 7 = 10^, 2 X 10^, 4 X 
10^, 10*, 2 X 10*, 4 X 10*. Time units are Larmor times for a 
7 = 10* particle . 



The asymptotic values of the difTusion coefficients, after 
the non-diffusive regimes, are plotted in figure |8] It is clear 
that K|| < K± < KBohm, in the energy range correspond- 
ing to fcniaxT'g ~ 1. In a uniform field, quasi-linear theory 
predicts kx < KBohm < However, in the highly non- 
uniform field produced by the non-resonant instability, the 
experimentally determined transport properties are much 
more complicated. Nevertheless, the parallel diffusion coeffi- 
cient, along with the shock speed, defines the residence and 
acceleration timescales. While the simulations presented in 
this paper do not include resonant wave excitation, which 
would act to reduce the diffusion coefficients even further, 
it is worthwhile speculating on the effects that our results 
would have on the shape of the CR spectrum and the rate 
of the sh ock acceleration pro c ess. T he acceleration rate de- 
fined by iLagage fc Cesarskvl (|l983l l is proportional to the 
magnetic field strength. Amplification of the magnetic field, 
as in the snapshot we selected, above interstellar medium 
values w ould therefore increa s e the maximum energy limit 
given by iLagage fc Cesarskvl (|l983l ). Inclusion of resonant 
wave excitation and application to the case of a circumstellar 
wind would increase that limit even further ()Bell fc Lucekl 
I2OOII ) as, of course, would a higher saturation level of the 
non-resonant waves. 

Turning to the shape of the spectrum, at a strong shock 
front with a compression ratio of four the standard result 
that / oc p~* applies provided that particles are isotropised 
and that their motion beco mes diffusive withi n a residence 
time either side of the shock (|Duffv et al.l[l"996[ l . In our sim- 
ulations the isotropisation is very rapid, on the timescale 
of one to two Larmor times. However, the parallel motion 
only becomes fully diffusive after almost ten Larmor times 
and is super-diffusive (< Ax^ >oc where P > 1) up to 
that point. It has been shown (|Kirk. Duffv fc Gallan^ll996l ') 
that if particle crossings happen on a timescale below that 
required for diffusive transport then the resulting spectrum 
will be a power law given by / oc p"'-^"''' provided the distri- 
bution remains close to isotropy. In our super-diffusive case, 
/3 > 1, this would imply a spectrum that is harder than the 
diffusive shock acceleration result. Therefore, our results are 




Figure 8. t||/KBohm (solid) and K,j_/KQahin (dashed) with 
KBohm evaluated in a 3/iG field, at t=6. 



consistent with the rapid acceleration of a relatively hard 
spectrum, but clearly more work needs to be done on res- 
onant wave excitation and overall dynamical range before 
drawing definitive conclusions. 



5 CONCLUSIONS 

We have reported on the results from high-resolution MHD 
simulations of the non-resonant current driven instability. 
The long term evolution of the field is still uncertain, as the 
non-linear development leads to very low density cavities 
in the absence of significant numerical dissipation. To over- 
come this, future simulations should include a non-linear 
feedback on the cosmic ray s perhap s using a hybrid code, 
similar to the work of iLucek fc Belli |2000. ) . or by perform- 
ing large-scale particle in cell ( PIC) simulations. To this end, 
iRiquelme fc Spitkovsk^ (12007] ) have performed 2D PIC sim- 
ulations o f the relativistic generalisatio n of the Bell-type in- 
stability (iReville. Kirk, fc Duffvl [20061 ) and preliminary re- 
sults show that the field saturates with both the perpen- 
dicular and parallel magnetic field well in excess of the ini - 
tial mean value. On the other hand, iNiemiec fc Pohli (|2007l ) 
have reported on 3D PIC simulations, with initial param- 
eters chosen such that the non-resonant mode is expected 
to be observed. However, the filamentation instability dom- 
inates the growth of the field and the perturbations do not 
grow much beyond the initial mean field value. The results of 
PIC simulatio ns of this instability remain inc onclusive. The 
simulations of IRiquelme fc Spitkovskvl ([20o3) are promising 
in that they do show the d evelopment of extended filaments, 
as predicted in lBelll (|2005l ). and this should be important for 
future investigations of particle transport. 

In recent years, magnetic field amplification has become 
a popular mechanism for accelerating cosmic rays beyond 
the Lagage-Cesarsky limit. We present here, the first at- 
tempt at investigating self-consistently the particle diffusion 
statistics in a self-excited magnetic field. The observed dif- 
fusion is anisotropic (k|| 7^ kx), with the low energy parti- 
cles' diffusion being largely determined by the field statistics. 
Thus at low energies > ft±. However, on length scales 
similar to the wavelength of the helical modes, motion par- 
allel to the shock normal is more diffusive. At high energies 
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the diffusion seems to be energy independent perpendicular 
to the initial field direction, but scales as 7^ parallel. This 
follows naturally since there are few waves with scalelength 
comparable to the gyroradii of these particles in our simu- 
lation box. However, the simulations do achieve sub-Bohm 
diffusion with only modest amplification of the effective am- 
bient field. 

It is cleax that in the presence of non-linear waves the 
diffusion properties of particles differ greatly from the stan- 
dard picture. Based on the results presented here it is dif- 
ficult to say if the super-diffusive or sub-diffusive regimes 
will have a significant effect on the shape of the spectrum. 
An estimation of the acceleration timescale of cosmic rays is 
beyond the scope of this paper, nevertheless we have demon- 
strated that self-excited turbulence does reduce the diffusion 
coefficient of the cosmic rays and is therefore likely to lead 
to acceleration beyond the Lagage-Cesarsky limit. 
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